      subroutine scmfdp1

      use param

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'xsec.h'
      include 'cmfdp1.h'
      include 'sanm.h'
      include 'sanmjin.h'
      include 'setls.h'
      include 'eigv.h'
      include 'nem.h'
      include 'scmfdp1.h'

      real(NBF) jnetl1,jnetr1,jnetl2,jnetr2

      call setcfls

      srcsrcd1=0
      srcsrcd2=0

      do iic=1,8
        if (mod(iic,2).eq.1) then
          iem=1
          srcsrc1=0
        elseif (mod(iic,2).eq.0) then
          iem=2
          srcsrc2=0  
        endif
        do im=1,ng
          do j=1,tnode
            tdl=dtl(iem,j,im)
            tdr=dtr(iem,j,im)
            hdl=dhl(iem,j,im)
            hdr=dhr(iem,j,im)
            czf(1,j,im)=-tdl+hdl
            czf(2,j,im)=-tdr-hdr
            if (iem.eq.1) then
              digf(j,im)=(tdr-hdr+tdl+hdl)+xsr(j,im)*hz(j)
            elseif (iem.eq.2) then
              digf(j,im)=(tdr-hdr+tdl+hdl)+(4*xsr(j,im)+5*xst(j,im))*de3*hz(j)
            endif
          enddo
        enddo

        do im=1,ng
          j=1
          del=digf(j,im)
          delinvf(j,im)=1/del
          jb=j
          do j=2,tnode
            deliauf(jb,im)=delinvf(jb,im)*czf(2,jb,im)
            del=digf(j,im)-czf(1,j,im)*deliauf(jb,im)
            delinvf(j,im)=1/del
            jb=j
          enddo
        enddo

! source term
        do im=1,ng
          do j=1,tnode
            ss=0
            do m=1,ng         
              ss=ss+avgflx(1,j,m)*xssm(j,im)%from(m)
            enddo
!            src(1,j,im)=reigv*xchi(j,im)*psif(j)+ss*hz(j)
            src(1,j,im)=reigv*xchi(j,im)*psif(j)+(ss+2*xsr(j,im)*psphi(2,j,im))*hz(j)
            srcsrc1=srcsrc1+src(1,j,im)*src(1,j,im)
            if (iem.eq.2) then
              src(2,j,im)=-2*de3*(src(1,j,im)-xsr(j,im)*psphi(1,j,im)*hz(j))
              srcsrc2=srcsrc2+src(2,j,im)*src(2,j,im)
            endif
          enddo
! forward substitution
          j=1   
          rhs=src(iem,j,im)
          y(j)=delinvf(j,im)*rhs
          jb=j
          do j=2,tnode
            rhs=src(iem,j,im)-czf(1,j,im)*y(jb)
            y(j)=delinvf(j,im)*rhs
            jb=j
          enddo            
! backward substitution
          j=tnode   
          psphi(iem,j,im)=y(j)
          jt=j 
          do j=tnode-1,1,-1
            psphi(iem,j,im)=y(j)-deliauf(j,im)*psphi(iem,jt,im)
            jt=j
          enddo
          do j=1,tnode
            if (iem.eq.1) then
              avgflx(1,j,im)=psphi(1,j,im)-2*avgflx(2,j,im)
            elseif (iem.eq.2) then
              avgflx(2,j,im)=psphi(2,j,im)
            endif
          enddo                  
        enddo ! for im    

        err1=abs(srcsrc1-srcsrcd1)/abs(srcsrc1)
        err2=abs(srcsrc2-srcsrcd2)/abs(srcsrc2)
        srcsrcd1=srcsrc1
        srcsrcd2=srcsrc2
        if (err1.lt.epsm3.and.err2.lt.epsm3) then
          print*,"iteration # of the scmfdp1=",iic
          exit
        endif
      enddo      

      call dumflx
      
! Determine the incoming partial current J_in
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          btm1=beta(1,jm1,im)
          bt1=beta(1,j,im)
          btp1=beta(1,jp1,im)
          btm2=beta(2,jm1,im)
          bt2=beta(2,j,im)
          btp2=beta(2,jp1,im)
          ppm1=psphi(1,jm1,im)
          pp1=psphi(1,j,im)
          ppp1=psphi(1,jp1,im)
          ppm2=psphi(2,jm1,im)
          pp2=psphi(2,j,im)
          ppp2=psphi(2,jp1,im)
!          sffdml=(btm1*avgflx(1,jm1,im)+bt1*avgflx(1,j,im))/(btm1+bt1)
!          sffdmr=(btp1*avgflx(1,jp1,im)+bt1*avgflx(1,j,im))/(btp1+bt1)
!          sfluxl(1)=sffdml+flcor(j,im)*(avgflx(1,jm1,im)+avgflx(1,j,im))
!          sfluxr(1)=sffdmr+flcor(jp1,im)*(avgflx(1,jp1,im)+avgflx(1,j,im))  
!          sfluxl(2)=sflux(2,j,im)
!          sfluxr(2)=sflux(2,jp1,im)
!          sfluxl(2)=0
!          sfluxr(2)=0
!          jin(1,j,im)=0.25*sfluxl(1)+0.5*jnet_l+0.3125*sfluxl(2)
!          jin(2,j,im)=0.25*sfluxr(1)-0.5*jnet_r+0.3125*sfluxr(2)
!          avgflx(2,j,im)=0.8*(4*jin(1,j,im)-2*jnet(1,j,im)-sflux(1,j,im))
! NEM P3 (using a pseudo flux)
          jnetl1=-dtl(1,j,im)*(pp1-ppm1)-dhl(1,j,im)*(pp1+ppm1)
          jnetr1=-dtr(1,j,im)*(ppp1-pp1)-dhr(1,j,im)*(ppp1+pp1)
          jnetl2=-dtl(2,j,im)*(pp2-ppm2)-dhl(2,j,im)*(pp2+ppm2)
          jnetr2=-dtr(2,j,im)*(ppp2-pp2)-dhr(2,j,im)*(ppp2+pp2)
          sffdml1=(btm1*ppm1+bt1*pp1)/(btm1+bt1)
          sffdmr1=(btp1*ppp1+bt1*pp1)/(btp1+bt1)
          sffdml2=(btm2*ppm2+bt2*pp2)/(btm2+bt2)
          sffdmr2=(btp2*ppp2+bt2*pp2)/(btp2+bt2)
          sfluxl(1)=sffdml1+flcf(1,j,im)*(ppm1+pp1)
          sfluxr(1)=sffdmr1+flcf(1,jp1,im)*(ppp1+pp1) 
          sfluxl(2)=sffdml2+flcf(2,j,im)*(ppm2+pp2)
          sfluxr(2)=sffdmr2+flcf(2,jp1,im)*(ppp2+pp2) 
!
          ajil(1,j,im)=0.25*sfluxl(1)+0.5*jnetl1-0.1875*sfluxl(2)
          ajir(1,j,im)=0.25*sfluxr(1)-0.5*jnetr1-0.1875*sfluxr(2)
          ajil(2,j,im)=0.4375*sfluxl(2)+0.5*jnetl2-0.0625*sfluxl(1)
          ajir(2,j,im)=0.4375*sfluxr(2)-0.5*jnetr2-0.0625*sfluxr(1)       
        enddo
        if (bcb.eq.1) jin(1,1,im)=0; ajil(1,1,im)=0; ajil(2,1,im)=0
        if (bcu.eq.1) jin(2,tnode,im)=0; ajir(1,tnode,im)=0; ajir(2,tnode,im)=0  
      enddo

      end subroutine
